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Abstract 

Details are presented of a recently announced exact solution of a model consisting of triangular 
trimers covering the triangular lattice. The solution involves a coordinate Bethe Ansatz with 
two kinds of particles. It is similar to that of the square-triangle random tiling model, due 
to M. Widom and P. A. Kalugin. The connection of the trimer model with related solvable 
models is discussed. 

1 Introduction 

The dimer problem is one of the classic models of statistical mechanics. A dimer in this context is a 
particle that occupies two neighbouring sites of a lattice. In the dimer-monomer model dimers and 
monomers (particles occupying one lattice site each) are placed on a lattice such that they cover 
all sites, without overlap. Equivalently the monomers can be viewed as empty sites; the lattice 
is then partly covered with dimers. This model was introduced to describe diatomic molecules 
adsorbed on a substrate Attempts have been made in vain to solve this model exactly, that is, 
to calculate its free energy. The special case that there are no empty sites (monomers) is known 
as the dimer problem. There the dimers cover the lattice completely and without overlap. This 
model has been solved for planar lattices, independently by Kasteleyn Q and by Temper ley and 
Fisher |^ . Their solution is based on the possibility to express the partition function of the model 
as a Pfaffian. For many planar lattices the dimer problem can also be solved by means of the 
Bethe Ansatz. On the honeycomb lattice for example it can be formulated as a five-vertex model. 
This is a special case of the six- vertex model, whose Bethe Ansatz solution is well-known [^[-pT]| . 
A review of the dimer problem is given in fl^ . 

Inspired by the solvability of the dimer model, we consider lattice coverings by trimers. A 
trimer is a particle that occupies three lattice sites. We only consider triangular trimers, which 
live naturally on the triangular lattice. As in the dimer model, we require that these particles 
cover the lattice completely and without overlap. Thus every lattice site is covered by precisely 
one trimer. Figure |l| shows a typical configuration. 



As will be shown in Subsection 2.2 this model has a structure of domains separated by domain 
walls. The domains are hexagonal, and the domain walls form a honeycomb network. Simi- 
lar domain wall structures are used to describe an incommensurate phase of a monolayer of a 
monoatomic gas adsorbed on a hexagonal substrate . The entropy of such a network is largely 
due to the "breathing" of the cells: it is possible to blow up a domain and simultaneously shrink 
its six neighbours, or vice versa. 

Hexagonal domain wall structures also occur in the square-triangle random tiling model [Q. 
For that model a coordinate Bethe Ansatz was found by Widom ||l5| . The resulting Bethe Ansatz 
equations were solved analytically in the thermodynamic limit by Kalugin | p6[ . An exact solution 
of the trimer model was announced in p^ ; in the present paper we describe its the derivation. 



Figure 1: A typical configuration of the trimer model. Each lattice site belongs to precisely one 
trimer. 



The solution is very similar to that for the square-triangle tiling, and we closely follow Kalugin's 
arguments. The outline of our calculation is as follows. A transfer matrix for the model is 
formulated. After the choice of a reference state two types of elementary excitations are found. 
They are closely related to the above-mentioned domain wall structure of the model. In order to 
diagonalise the transfer matrix a coordinate Bethe Ansatz is set up in terms of the elementary 
excitations. The resulting semi-grand canonical ensemble is discussed. In the thermodynamic 
limit the Bethe Ansatz leads to a set of two coupled integral equations. These can be solved in 
a special case. From their solution the relevan t ph ysical quantities are computed. The results of 



the calculation are summarised in Subsection 3.5. We then consider the entropy as function of 



the density of down trimers. The model undergoes two phase transitions in the density of down 
trimers. 

Finally we discuss the relation between the trimer model, the square-triangle random tiling 
model, and yet another solvable model with a hexagonal domain wall structure. 



2 Preliminaries 
2.1 Sub-lattices 

Figure ^ shows a very regular configuration of the model, in which the trimers are positioned on 
a sub-lattice of the triangular faces. There are six such sub-lattices, which we number 0, 1, ... , 
5 as indicated in the figure. Note that the even-numbered sub-lattices consist of the up triangles 
while the down triangles constitute the odd-numbered sub-lattices. For a given configuration let 
A'^ denote the total number of trimers and let Ni denote the number of trimers on sub-lattice i. 
We wish to compute the entropy per trimer as a function of the sub-lattice densities 

These densities satisfy the obvious linear constraint 

PO + Pi + P2 + P3 + P4 + P5 = 1- (1) 



In Subsection 2.5 it will be shown that when toroidal boundary conditions are imposed the densities 



also satisfy a quadratic constraint: 

PaP2 + P2Pi + PiPo = P1P5+ P3P5 + PbPi ■ (2) 

Hence of the six sub-lattice densities only four are independent. In order to be able to set up a 
transfer matrix we pass to the grand canonical ensemble. The trimers on each sub-lattice i are 
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given a fugacity Wi = exp(/ii). After the transfer matrix has been diagonalised we shaU Legendre 
transform back to the canonical ensemble. 




2.2 Domains and walls 

Occupying sub-lattice completely while leaving the other five sub-lattices empty results in the 
configuration of the model shown in Figure ^. This arrangement does not admit local changes. 
However, it is possible to fiip a whole line of trimers. Such a line can be viewed as a wall 
separating two domains consisting of trimers on sub-lattice 0. These domain walls come in three 
types (orientations), corresponding to the three odd- numbered sub-lattices. When two walls of 
different type meet a wall of the third type is formed. A trimer on sub-lattice 2 or 4 occurs 
when three domain walls of different type meet in a Y, but this does not happen at an upside- 
down Y. Figure ^ show examples of how the three types of domain walls can meet. In a general 
configuration the domain walls form a hexagonal network. 

2.3 Transfer matrix 

In an allowed configuration of the model each lattice site belongs to precisely one trimer. This 
trimer sits either on one of the three lattice faces above the site or on one of the three faces below 
the site. Label the site with a "spin" t or | accordingly. 

Consider a horizontal row of lattice sites and assume that the trimer configuration below that 
row is given. It determines the spins on that row. The sites occupied by a trimer below have a 
spin I, while those not occupied by such a trimer must needs carry a spin "f. Now consider the 
next layer of lattice faces, above this row. In order to decide what trimer configurations on this 
layer are possible, it is sufficient to know which sites are already covered. This is precisely the 
information encoded by the spins. 

This shows that the model can be described in terms of a transfer matrix that connects two 
consecutive rows of spins. Let cr denote be the spin configuration on the lower row and t the spin 
configurations on the upper row. Consider all trimer arrangements (without overlaps) on the layer 
in between that are compatible with the spin configurations cr and t. (Generally there is at most 
one such arrangement.) The sum of their Boltzmann weights is the transfer matrix element T^-a-. 

The rows of lattice sites in the model come in two types, call them A and B, that are shifted 
with respect to each other. The rows of type A and B alternate. Hence there are in fact two 
transfer matrices. Tab for layers with upper row of type A and lower row of type B, and Tba for 
layers with upper row of type B and lower row of type A. The products TabTba and Tba^ab 
are double transfer matrices that act between rows of equal type. Take a lattice consisting of 2M 
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Figure 3: The configuration from Figure ^ admits line excitations. These domain wahs can meet 
in Ys (top) and upside-down Ys (bottom). The Ys are chiral; the mirror image of the Y shown 
here contains a trimcr on sub-lattice 4 instead of 2. The upside-down Ys are achiral. To guide the 
eye the trimcrs not on sub-lattice are shaded lighter; the numbers indicate their sub-lattices. 



layers and impose periodic boundary conditions in the vertical direction, that is, identify the lower 
and upper row. As usual the partition function of the model on this lattice is 

Z = Tt (TabTba)*' - Tr (TbaTab)*'. 

In the limit that M tends to infinity the partition sum is dominated by the largest eigenvalue of 
TabTba or TbaTab- 

The matrices Tab and Tba can be combined into a single matrix 



T : 

acting on vectors 



Tab 
Tba 



V'A 
V'B 



where ipA and ipB are "wave functions" on the rows of type A and B, respectively. If such a vector 
is an eigenvector of T with eigenvalue A then ipA and ?Ab are eigenvectors of TabTba and TbaTab, 
respectively, with eigenvalue A^. 

Some identification of the row types A and B could have been chosen in order to avoid the 
complication discussed. This amounts to skewing the lattice or, equivalently, the transfer matrix 
direction. Because it breaks a mirror symmetry of the system we have avoided this solution. 

2.4 Conserved quantities and elementary excitations 

In the configuration obtained by fully occupying sub-lattice 0, each row of spins consists of re- 
peating blocks tit- Therefore we group the sites into blocks of three, as in Figure |[ Number the 
blocks in a row from the left to the right. 
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Figure 4: Spins for the configuration from Figure For clarity the edges of the triangular lattice 
have been largely omitted. 



Consider a trimer configuration on a layer of the lattice. Let L denote the number of blocks 
per row and let no, . . . , denote the number of trimers on each sub-lattice. The horizontal 
and vertical lattice direction are viewed as "space" and "time" , respectively; the lower and upper 
row of the layer then are time-slices at times t and t + 1. Counting the number of t spins in the 
lower row, distinguishing by the position inside the block, gives: 



na+ni+ 712, 



(*) 

= "-2 + "3 + "4, 
(*) I 1 



From this we get 



(t) 
(t) 



n 



it) 
(t) 



L — no - ni + + 714, 
,1 ".T. — L - no + n2 + ns - 125. 
Counting the number of J, spins in the upper row gives: 



(t+i) 
i" 

(t+i) 
l» 

(t+i) 
I 



n3 + n4 + 715, 
n5 + no+ 711, 
ni + 77-2 + 713. 



From this we get 



(t+i) , (t+i) r 



i -r n;^, ' = L-no 



"2 



"3 
"3 



- 774, 
77,5. 



(3) 
(4) 



(5) 
(6) 



Comparing with (||) and (Q) with we see that the quantities 77|,, +77,|, and n„ i +77,|, 

are conserved between rows. 

These conserved quantities are nonnegative. The only row of spins for which both are zero 
consists entirely of blocks tit- There is only one way to fit a layer of trimers above this row. Of 
course the row of spins above that layer consists again entirely of blocks tit- This row state will 
be chosen as the "empty" or reference state for the Bethe Ansatz in Section ^j. 
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A row of spins with rij,, +nm]» = 1 and n„j = is obtained by replacing one block, say 

at position x, in the reference state with J,||. There is only one possible configuration of trimers 
on the layer above, see Figure |. The row above consists of blocks "f If except for one block lit at 
position a; — i. Thus the transfer matrix has shifted the block iit in the lower row half a step to 
the left in the upper row. This block is a left-moving elementary excitation of the reference state. 
It will be called an L-particle. Similarly the block ||| is an elementary right-moving excitation, or 
R-particle. The conserved quantities n|,, +n,i, and n„| are the number riL of L-particles 

and the number tir of R-particles, respectively. 









J) ( 























Figure 5: There is only one way to fit trimers below a row consisting of one block amidst 
blocks m. It leads to another such row of spins below. 



The particle content of the blocks m, m and m has now been determined. For each of the 
other five blocks both til and riR are greater than zero. Therefore these blocks are combinations 
of the elementary excitations. They will be discussed in more detail in Subsections 3.3-3.5. 

We have found no other conserved quantities than riL and riR (except in the special case when 
riL = or npi — 0). 



2.5 World lines and quadratic constraint 

Divide the lattice into hexagonal patches containing one face from each sub-lattice, in such a 
way that the lower middle triangle of each patch belongs to sub-lattice 0. There are six trimer 
configurations possible on such a patch. Decorate each patch according to this configuration as 
shown in Figure ^ It is tedious but straightforward to verify that the decorations of the patches 
making up the lattice fit together, giving a set of solid and dashed lines running from the bottom 
to the top of the lattice. It can also be checked that the crossings of these lines with the lattice 
rows correspond to the locations of the L-particles (solid lines) and R-particles (dashed lines). 
Hence these lines are the "world lines" of the L-particles and R-particles where the horizontal and 
vertical lattice direction are viewed as "space" and "time" respectively. 




Figure 6: Decompose the triangular lattice into hexagonal patches such that the lower middle 
triangle of each patch belongs to sub-lattice 0. The other triangles, in counterclockwise order, 
then belong to sub-lattice 1, 4, 3, 2 and 5. These patches can be decorated with the world lines 
of the L-particles (solid) and R-particles (dashed). 

Impose toroidal boundary conditions. We now derive the quadratic constraint (|^) by the same 
method as that used for rectangle-triangle random tiling models in fl^ and ||l^. Cut the torus 
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open along a horizontal row of sites, so that the model is now on a cylinder. By stacking a number 
of copies of the configuration on top of each other we can achieve that each world line winds 
around the cylinder an integer number of times. Let 2M be the number of rows. In each row we 
can count the number of L-particles and R-particles: 

TiL = i — no - 111 + ?i3 + 714, 
riR = L - no + ^2 + ^3 - ns. 

Summing over the entire lattice yields 

2MnL = 2LM ~ - Ni + N^i + N^, (7) 
2MnYi,^2LM-No + N2 + N:i~N^. (8) 

Since the world lines of the L-particles do not cross each other, they all have the same winding 
number Wl- In order to compute this winding number consider the total leftward movement of 
the L-particles. There are nL such particles, each of them winds around the cylinder Wl times, 
and each winding constitutes a movement of L blocks to the left, so the total leftward movement 
amounts ni^Wi^L blocks. It can be seen from Figure ^ that an L-particle moves half a block to 
the left at a trimer on sub-lattice 2 or 5, while there is no horizontal movement of L-particles at a 
trimer on another sub-lattice. Summing over the entire lattice shows the total leftward movement 
of the L-particles to be ^{N2 + N^) blocks. Hence 

nLWLi = ^(A^a + iVs). (9) 

Fully analogously one has 

nnW^L = ^{Ni + N^) . (10) 

Now consider the number of crossings of L-particle lines and R-particle lines. There are nL L- 
particles each winding leftward around the cylinder Wl times, and riR R-particles each winding 
rightward around the cylinder Wr times, so the number of crossings is nLnR(WL -t- Wr). From 
Figure ^ it is seen that crossings occur precisely at trimers on sub-lattice 2 or 4, so that the number 
of crossings is N2 + N4. Equating these two expressions for the number of crossings gives 

nLnR(WL -I- Wr) =N2+ N4. 

Substituting into the above equation first (^) and then (0) and (|), then muhiplying by 2LM 
and finally using 

2LM = Nq + Ni + N2 + N3 + Ni + N5 

yields 

N1N3 + + N5N1 = NqN2 + N2N4, + N^Nq. 

Dividing by N"^ gives (|^) . For the sake of the argument we have stacked a number of copies of the 
original configuration on top of each other, but the resulting (||) is not affected by this. 

3 Bethe Ansatz 

In this section we describe a Bethe Ansatz (BA) that diagonalises the transfer matrix. Since the 
particle numbers nL and nR are conserved quantities the transfer matrix is block diagonal. We 
begin by considering the sector with nL = and nR — and then pass to sectors with higher 
particle numbers. 
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3.1 No particles 



The only state in the sector til = 0, tir = is the reference state which consists entirely of 
blocks tit, so this sector is one-dimensional. Therefore the transfer matrix acting on this sector 
is trivially diagonal. The layer between two consecutive rows in the reference state consists of L 
trimers on sub-lattice 0, so its Boltzmann weight is . It is the eigenvalue of the transfer matrix 
in this sector. For convenience we define the "reduced" transfer matrix T to be the transfer matrix 
T divided by . 



3.2 One L-particle 

Consider a row of spins containing a single L-particle (lit) at position x. The transfer matrix has 
shifted this particle from position x -I- ^ in the row below half a step to the left, see Figure I The 
layer between the two rows contains L—1 trimers on sub-lattice and one trimer on sub-lattice 5. 
Hence the action of the (reduced) transfer matrix on the "wave function" is given by 

{fi^){il]x) = ^4,{iiU + \)- 

Wo 

(We use ttt X as notation for the row configuration that has a block lit at position x and blocks 
tit at the other positions.) The solution of the eigenvalue problem Tip = Aip is 

iP ill^ x) ^ A^u"", 

where is some constant, and 

A = — u2 . 

Wo 



3.3 One L-particle and one R-particle 

Consider a row of spins containing an L-particle (ttt) at position x and an R-particlc (tii) at 
position y, with x < y. If the particles are apart, this situation was formed by shifting the 
L-particle to the left and the R-particle to the right: 

{maiU,niy)^^^i^{iiU + hniy~^) ity-x>i. (n) 

Wq z z 

(We write the arguments of ip in order of increasing position; for example, the notation in the 
LHS of ( pT| ) implies that x < y.) If however the particles are next to each other, the situation was 
formed from a "bound state" (iii), see Figure || (top): 

[mm z - i tii ^ + i) - ™ vxiii z). (12) 

This bound state was formed from another type of bound state (ttt)- 

{fmu z) = ^^vxttt z-\)-r ^^vxttt z + i). (13) 

Wo ^ Wo 

The two terms correspond to two chiral configurations one of which is depicted in Figure |. This 
bound state was formed from an R-particle and an L-particle next to each other, the R-particle 
sitting to the left of the L-particle: 

(f V')(ttT z) = — VXtii Z-I UU + h)- (14) 

Wo 

This configuration can have arisen from the same bound state again. The alternation of this 
bound state and the situation where the R-particle and L-particle are next to each other (tii iit) 
corresponds to the vertical domain wall in Figure ^. The configuration where the particles are 
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next to each other can also have arisen by shifting the R-particle half a step to the right and the 
L-particle half a step to the left, see Figure || (bottom): 



(f V')(Tii z-\,ii]z + \) = w^i^im z) + ™^(Tii z-l,ii]z + 1). (15) 



Finally, a configuration where the particles are apart was formed by shifting the R-particle half a 
step to the right and the L-particle half a step to the left: 

{fmily,llU)^^^^{niy-\AlU + \) ifx-y>i. (16) 

We want to solve the eigenvalue equation Tij) — A^/; for (pl])-(p^. The eigenvalue equation for 
(|ll|) and (|l|) is satisfied by 

^/'(iii z) = AuvU^'v^ 



with eigenvalue 



r I wi _i 

A = — M2 V 2 . 

Wo wq 



Similarly the eigenvalue equation for ( |14[ ), (|15|), and (|16D, with the same value for A, is satisfied 
by 

V'(Tii J/,iiT a:) = A„„w^u^ \ix-y>l, 



where 



tU v] , (17) 



2„,,2 



uv 



D = ^. (18) 
W1W5 

The eigenvalue equation for ( p^ is satisfied too if 

A 

where 

Suv = — u -I V 1 V . (19) 

W1W5 \W5 Wi J \ wjwf / 

The above analysis suggests to interpret the bound state iii as LR (in that order) and the bound 
state TTT as RL. The eigenfunction is then written 



ipiL x,Ry) = AuvU^'vy, 

■0(R ?/, L x) = B^uA^uvyw^ ifa;-y = l, 
I DByuAyuV^W" if x-y = 0. 
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3.4 Two L-particles and one R-particle 

A similar but more tedious analysis can be carried out for the sector with two L-particles and one 
R-particle. There is a new bound state (itt) that can be interpreted as LRL. A solution of the 
eigenvalue problem Tip = Atp is given by 



ipiL xi,R y,L X2) = Bvu-...A 



V(2) 



Xi V X2 



if 2:2 - y > 2, 
if 2^2 - y = 1, 



V'(R 2/,L xi,h X2) 



Stt ^''"ir(l)'",r(2)'''^'",r(l)'"7r(2) Xi — y >2, 

X^TT -^"«x(l)^f"x(l)«,r(2)^^"7r(l)"7r(2) if^^l— 2/=!, 
.X^TT ^-^W«x(l)^fMx(l)M,r(2)^^"7r(l)"7r(2) if ^^l ^ J/ = 0, 

where tt runs through the permutations of {1,2}. The amplitudes must satisfy 



= -1 



A 


A 


A 


A 


A 


A 

-^U^f UiV 







(20) 
(21) 



with Suiv given by dl^). Note that the amplitude ratio in ( pO| ) does not depend on v and that the 
amplitude ratio in ([2l[) does not depend on Uii . The eigenvalue is given by 



T W5 i W5 i Wi _i 

A = — ul — m| — V 2 . 

Wo Wo Wo 



3.5 Arbitrary particle numbers 

With two L-particles and two R-particles, there is a new bound state (iti) that can be interpreted 
as LRLR. This completes the list of possible blocks and their interpretation in terms of particles, 
see Table |[ 



Table 1: The three-spin blocks. 



spins 


particles 


Tit 


none 


UT 


L 


Ui 


R 




LR 


TTT 


RL 


itt 


LRL 


tn 


RLR 


in 


LRLR 



The solution given above of the eigenvalue problem — for two L-particles and one 
R-particle generalises to the higher sectors. Before describing this generalisation we introduce 
a notational convention. The index i, running from 1 to til, will be used to number L-particle 
positions and Bethe Ansatz. The index j, between 1 and riR, will refer to R-particles. Now 
consider a succession of L-particles with coordinates xi < X2 < ■ ■ ■ < Xm^ and R-particles with 
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coordinates yi < 2/2 < ■ ■ • < J/rin ■ (Note that Xi — x^+i can arise only from a block LRL or LRLR, 
so Xi — yj = a^i+i for some yj.) The value of the wave function is given by 

Vjlparticle sequence) = X! X! 11 ^'^'^ ^- -^ ^ ■■ IT ("^(0 ^'^'^ ""^0)) ' ^'^'^^ 

TV a 

where tt and a run through all permutations of {1, 2, . . . , til} and {1, 2, . . . , tir}, respectively. We 
shall now describe the factors in the RHS of (p^). For each segment R , L Xi in the particle 
sequence with Xi — yj ^ 1 there is a factor By^i^.^u^^iy For each such segment with Xi — yj = 
there is a factor DB^^^.^^^^.y The amplitude A... depends on the sequence of the variables u 
and V corresponding to the sequence of L-particles and R-particles. The u are in the order 
'"7r(i), Uir(2)i ■ ■ • iW7r(nL) ^'^d the V are in the order Va(i)TV„^2), ■ ■ ■ ji'cr(nri), but the two sequences 
interlace. The amplitudes A... are defined up to an overall factor by the conditions 



A 



A 

■^...U^lUi... 



A 



A 



— SuiVj 



with SuiVj given by (|l9|). Finally comes the product of all the and w^q-j. The eigenvalue for 
the eigenfunction is given by 

A=n^4n^.-i (23) 

1=1 J=l 

We have no rigorous proof that the above solution is correct for all sectors, but using computer 
algebra we have verified it for t^l + < 5. 

It should be noted that the formulation of the solution depends on the particle interpretation 
of the three-spin blocks. The particle content of each three-spin block is determined by til = 
film* + JT-.f* and npi — -I- ri,|,, but the order of the particles within a block can be chosen. 
For example, we could interpret J,tt as LLR, LRL or RLL. The choices in Table |l| lead to a 
simple description of the eigenfunctions; each factor D or B depends only on two successive 
particles. Other choices than those in Table Q would make the formulation of the eigenfunctions 
more awkward; there would be more factors than just D and B, and some would depend on 
non-successive particles. 

3.6 Bethe Ansatz equations 

Impose periodic boundary conditions in the horizontal direction. This means that the wave func- 
tion must not change if the leftmost particle (L at xi or R at j/i) is moved to the corresponding 
position at the other (right) side of the system: 

V'(L xi, . . .) ^ . . ,L xi + L), 
^{R yi,...) =7A(...,Ryi-f L). 



The eigenfunction ip given by (22) satisfies these conditions if the Bethe Ansatz equations (BAEs) 
hold: 



uf = i~r--'Y[S^^,^, (24) 

^L^(_)nH-l-Q5-^l^. (25) 
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Note that although the description of an eigenfunction in terms of the u and v involves factors 
(|l^) and (|l|), the BAEs only contain factors 
Upon substitution of 



\w1W2wl ) 



the expression ( Jl^ ) for S^v simplifies to 



and 



(26) 



The BAEs (E 



and (p^) can then be written 



1 + 



/wgZ«3W4 Y / W1W3W5 
\W1W2wlJ \W0W2W4, 

I W1W3W5 
WQW2W3 J \W0W2W4 



3 

'W1W4W5 \ 



- Vj 

\L+riR-l TT C-1 ^ 



(27) 



(28) 



These equations can be considered the key result in the exact solution of the model. They deter- 
mine the possible values for the ^ and the the rj. These in turn determine the eigenvalues and the 
eigenfunctions of the transfer matrix, 



A = 



\W0W1W2 J \W0W4W5 J 



z=l 3 = 1 



(29) 



where we have reintroduced the factor Wq that was omitted as of Subsection 3.1 



As a check on the Bethe Ansatz we determined the eigenvalues of the transfer matrix for small 
system size by (brute force) numerical diagonalisation; the same eigenvalues were obtained by 
numerically solving the BAEs. 



4 Thermodynamics 

We are interested in the behaviour of the model as a function of the sub-lattice densities, that 
is, the canonical ensemble. In order to set up a transfer matrix we have passed to the grand 
canonical ensemble, which is controlled by sub-lattice weights (or chemical potentials) instead of 
sub-lattice densities. In this section it turns out that the transfer matrix leads to a semi-grand 
canonical ensemble. It is controlled partly by densities (essentially the two conserved quantities) 
and partly by chemical potentials. We describe the Legendre transformation from this ensemble 
back to the canonical ensemble. We also look into the symmetries between the sub-lattices and 
how they appear in the semi-grand canonical ensemble. 

4.1 Legendre transformation 

In passing to the grand canonical ensemble each trimer on a sub-lattice i was given a weight 
Wi — exp(/ii). Certain combinations of these weights occur in the BAEs (^7|) and ( psj ) and in the 
expression (^9|) of the transfer matrix eigenvalue in terms of the BA roots. It is convenient to 
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assign names to the corresponding combinations of the chemical potentials /i^ = log Wi, 

<^R = ^ [(3^0 - 3pi + /i2 + M3 - M4 - Ms) + PL (-Mo + Ml - M2 + M3 - M4 + Ms)] , 
ML = ;^ (-Mo - Ml - M2 + M3 + M4 + Ms) , 
MR = ^ (-Mo + Ml + M2 + M3 - M4 - Ms) , 

where 

ML = and pR = — 

denote the densities of the particles L and R. With these definitions the BAEs ( p7| ) and ( p8| ) can 
be written 

(e^^e.)^ = (-)-+--^ n (30) 

V + Vi 



i=i 



(e'-V^)' - (-)^+"^-^nC^^^, (31) 

r=i + 



while the eigenvalue expression (^9|) becomes 



/ "L "R \ 2 

A = exp(L/io + n-LML + ?^rMr) !!("'':') ) • (3^) 

Taking the logarithm, dividing by L, and letting L to infinity gives the free energy per trimer in 
the thermodynamic limit: 

^^(ml, Mr; Mo, Mi, ■ • ■ , Ms) = *(ml, Mr; 0l, 0r) - MlMl - MrMr - Mo, 

where 

^ / "L "R \ i 

$(pl,Mr;0l,</'r) --^lini^-^logl n(-'^j)) ' (^3) 

It is the free energy in a semi-grand canonical ensemble where the numbers of trimers on the 
different sub-lattices may vary but are subject to the constraints imposed by fixing the particle 
densities 

Ml = 1 - Mo - Ml + M3 + M4, (34) 
MR = 1 - Mo + M2 + M3 - Ms- (35) 
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In order to do the Legendre transform to the canonical ensemble the derivatives of Q with respect 
to fiQ, ^i, fi^ have to be taken. This gives the ensemble average densities 



Po = 



Pi 



P2 



P3 



Pi 



P5 = 



PR 



PR 



PR 



2 " 



^PR 



^PR 



9$ 



PL 



PL 



PL 



2PL 



^PL 



d(j)R 

90R 



1 

4 
1 
4 
1 
4 
1 

4^^ 



PL 



PL 



PL 



1 

4 
1 
4 
1 
4 
1 

4PR, 



PR 



PR, 



PR, 



4PL - ^PR, 



4PL - ^PR. 



1. 



(36) 
(37) 
(38) 
(39) 
(40) 
(41) 



In Subsection 2.1 it was seen that because the sub-lattice densities satisfy two constraints, four of 
them are independent. Equations (^^-(^Tj) express the sub-lattice densities in terms of only four 
quantities, namely pL, PR, and Therefore these four quantities must be independent, 

and the sub-lattice densities given by (^-(^ must satisfy the two constraints, (|l|) and (||). This 
can also be verified by direct computation. The entropy per trimer is 



S{po,pi, . . . ,P5) = -0+ Vp£/i£ = -$ 



9$ 



/ , I 

=0 



5R- 



(42) 



It is remarkable that the chemical potentials /iq, pl and /ir that occur in the expression ( p^ ) for 
the eigenvalue have disappeared in the Legendre transformation. As a consequence $ and hence 
the densities po, Pi, • • ■ , Ps and the entropy S are now functions of four parameters: the particle 
densities pL and /jr and the potential- like quantities (jfh and 0r. These are just the parameters 
that govern the BAEs (H) and (|l|). This agrees with the fact that the canonical ensemble also 
has four parameters. 



4.2 Symmetries of the parameter space 

For the reference state of the BA sub-lattice was chosen. Since the model is invariant under 
horizontal translations over a single lattice edge sub-lattice 2 (or 4) could have been chosen instead. 
The original situation can be regained by renumbering the sub-lattices i i — 2{ mod 6). The 
sub-lattice densities p[ in the new numbering are related to the densities pi in the old numbering 
by 

Po = P2, P'l = P3, etc. (43) 
and analogously for the chemical potentials pi. From this one computes 

Pl = 2-Pr, Pr = 1 + /?l-Pr, = -(/"l + </'r, = 

Similarly the model is invariant under reflection in a horizontal line. The corresponding sub-lattice 
renumbering is i ^ i + 3( mod 6). This gives: 

Pl = 2 - PL, Pr = 2 - PR, = (/iL, (t'n = 0R- 

The model is also invariant under reflection in a vertical line. For the line passing through sub- 
lattices and 3 the renumbering is i —^ —i{ mod 6). Obviously this is nothing but interchanging 
left and right, so 

Pl = PR, Pr = Pl, = 0r, = 



14 



Together these three transformations generate a group of order twelve. In Subsection 6.3 we 
shall find four "families" of points in the parameter space where the entropy of the model can 
be computed exactly. These four families turn out to be related by symmetries from this group. 
Note that under this symmetry group the free energy Q and the entropy S are invariant, so $ 
transforms in a certain way. For example, for the translation (^3|) the transformation is 

^■(2 - PR, 1 + PL - PR,; -0L + 0R, = *(PL, PR-, 0L, (/-r) " ^^L. 

Finally the model is invariant under some rotations. As an example, consider the rotation over 
2tt/3 about an up triangle of the lattice. The sub-lattice renumbering is: 1 ^ 3 ^ 5. This does 
not give a simple transformation of /Ol, Pr, (f'L and (j}R. This can be explained as follows. In the 
definition of these four parameters the direction in which the transfer matrix acts plays a special 
role. Rotations do not preserve this direction, in contrast to the translation and the two reflections 
described above. The symmetry group generated by all these operations is of order 36. 



5 Integral equations 

In Section ^ two sets of BAEs were derived. These equations can be solved numerically, for system 
size L up to a few hundred, say. This can be done essentially in the full parameter space. (The 
regions where numerical complications arise can be mapped to regions without such difficulties 



by means of symmetries from Subsection 4.2.) We however want to get analytic expressions for 
the physical quantities of the model in thermodynamic limit. In the present section the BAEs 
in the thermodynamic limit are turned into two integral equations for two complex functions. 
These functions are multivalued, and their monodromy properties are obtained from the integral 
equations. The functions are then determined from their monodromy and analyticity properties. 
In the next section these functions will be used to compute physical quantities of the model. 

5.1 Derivation 

We shall now in the usual fashion derive integral equations from the BAEs (|3^) and (^l|). The 
logarithmic version of these BAEs is 

LFi^i^,) = (riL + riR - l)7ri (mod 27ri), (44) 
LFr{7Jj) = {L + nR- l)m (mod 2m), (45) 

where 



Fl (z) ^^ogz-jJ2 - ^j) - log(^ + ^7')] + + ^ log rj, , (46) 

Fr,(z) =logz--J2 [log(^ - e.) - log(z + C')] + <^R + 7 E (47) 



i=l i=l 



The derivatives of these functions are 



/LW = --7Ef^ -^), (48) 

For the understanding of the structure of the solutions to the BAEs we rely on numeric compu- 
tations for finite system size. For many values of the parameters pL, PR, '/'l and 4>r the BA roots 
for the largest eigenvalue show the following features. The roots and -qj lie on smooth curves 
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in the complex plane. When the system size becomes large these curves tend to well-defined limit 
shapes. These limit curves will be called S and H. The sets {^^j and {rjj} (and hence also the 
curves S and H) are invariant under complex conjugation; this implies that 



and 



/r(^*) - /r(^)' 



(50) 



The curve 5 crosses the positive real axis, whereas H crosses the negative real axis. Note that by 
( p^ ) the roots are solutions of 



LFl(0 = ("L +nn- l)7ri (mod 27ri). 



(51) 



This equation defines discrete points on the curve Rei^L(0 = 0. The roots occupy a succession 
of these points, without holes: 

i[i^L(6+i)-i^L(C»)] = 27ri. (52) 

By holes we mean solutions of ( ^ lying between and curve Rei^L(0 ~ that are 

not contained in the set Similary for {r/j}: 



L[Fn{vj+i)-Fn{vi)]^27ri. 



(53) 



Let 6l and 6r denote the end points in the upper half plane of S and H. When the system becomes 
large and tend to and 6l, respectively, while 771 and ?7„j^ tend to &r and b^, respectively. 
Figure]^ shows the distribution of the roots for the largest eigenvalue in a given sector riL, n^. 



2i 



•-2 



-1 



-2i- 



. 2 



Figure 7: Distribution of the BAE roots for the largest eigenvalue. The ^ are on the right, the rj 
on the left. The parameters have the values (pi^ = —0.46, (pn = —0.653, nL — 15, = 18 and 
L = 30. 

We assume that the condition that there are no holes in the sets of roots {^i} and {r/j} also 
holds in the thermodynamic limit. There (^2|) and (|5^ ) can be written 



LfRiVj)iVj+i - Vj) = 27ri, 



so that the sums in (48) and (^9|) can be turned into integrals 



/r('7) d??, 



^h(^)4 + 2^I(c^ + FTTi)^^(^)'^^- 



(54) 
(55) 



(56) 
(57) 
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From these equations it is seen that /l(z) has branch cuts H and — and that fn{z) has branch 
cuts S and — From the same equations it is easily computed that zfi^{z) and z/r(z) are 
invariant under z ^ —z^^. Therefore we substitute 

z-z-^^z, (58) 

and define 

zh{z)^gi,{z) and zfj^{z) = g^{z). (59) 

The two branch cuts H and — H^^ of /l(^;) then collapse to a single branch cut H of 5l(^), and 
similarly for /r(2:). The equations (|56| ) and (^^ become 

3l(^) = ^ + 7^ f -A-r gnifl) dry, (60) 
27ri JiiTi- z 

^^^^^ = ^ + i I |~ -"^^^'^ "^^^ ^^^^ 



The functions fhiz) and /r(z) and hence also 5l('Z) and 5R(i) contain all information about the 
BA roots and rjj that is needed to compute the densities pL and pR, the phases 0l and c/fiR, and 
the semi-grand canonical free energy $. 

The integral equations (^0|) and ( |6l| ) are very similar to the equations obtained by Kalugin jl^ 
for the square-triangle random tiling model. He tackles his equations by exploiting the monodromy 
properties of the functions. We shall use the same method for our integral equations, closely 
following Kalugin's argument. 

5.2 Monodromy properties 

The first step in solving the integral equations (^^ and (^l|) (but only for a special case to be 
defined below) is to establish the monodromy properties of the functions 5l(^) and gniz). From 
( |60| ) it is seen that the contour H is a branch cut of the function gh{z). Consider this function 
on a path Fjj that crosses the contour H in some point zo, as in Figure ^(a). The jump of the 
function gi^ over the contour is 

5L(io)after " 5L (2o)before = TT^T </ gRiv) d?? = ffR(zo)- 



27ri 



|=£ V- Zo 



Hence the analytic continuation of gi^{z) along the path Fjj through the contour H is gi^{z) — gjx.{z) . 
From ( |6l| ) the contour H is not a branch cut of the function g^{z), so the analytic continuation 
of gniz) along Fg is just gniz). Therefore the effect of analytic continuation along Fjj on a linear 
combination ai^gi,{z) + a^^g^{z) is given by 



It can be derived analogously that the monodromy operator for the path T~ in Figure |^(a) is given 
by 

The operators Fjj and T~ generate the full group SL(2, Z). Now consider the special case that the 
end points and 6£ of 5 coincide with the end points ^r and b'^ of H, and that the contours do 
not meet in other points. Then the paths Fjj and T~ are no longer defined, but their composite 
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r^jFe is, see Figure ^(b). Since this operator is of order six ghiz) and (7r(z) are single- valued 
functions of the variable 



z-b 
z ~ h* 



,1/6 



(62) 



where 6 = 6l = &r and 6* = = 6^ are the common end points of the contours. The inverse 
transformation is 

b*t^ - b 
z — — ^ . 

<6 - 1 

Since (rjjre)gR(i) = 5l(^), the functions ghiz) and gn{z) are different branches of a single 
function g{z). 





(a) 



(b) 



Figure 8: The contours ^ and H. (0l 
PL = 0.76, PR = 0.93.) 



-0.46, 



-0.653 and (a) pL = 0.5, pr = 0.6, (b) 



In the remainder of this paper we shall, unless stated otherwise, consider the case that the 
end points of the contours 5 and H coincide, and that these are the only common points of these 
contours. 



5.3 Analyticity properties 

It follows from (^^ that gh{z) is analytic everywhere except on the branch cut H. Similarly gB.{z) 
is analytic everywhere except on S. In particular it is analytic on the contour H, except perhaps 
at the end points, as these lie also on S. It then follows from ( |60| ) that ghiz) remains finite if 
z approaches a point (not an end point) on its branch cut H. An analogous statement holds 
for gB.{z). To summarise, ghiz) and gniz) are finite everywhere except perhaps at b and b*. 

It was derived above that they are branches of one function g{z), which in turn is a single- 
valued function of t. Fix this function h{t) by choosing that at i = e'^'/'^ it corresponds to gh{z) 
(at z = oo). Figure ^ shows where in the t-plane all the branches of g(z) = h{t) are situated. Note 
that the branch at t = e~^^/^ is 5r(z) (at z — co). 

Since g{z) is finite everywhere except perhaps at z = b and z = b*, h{t) is analytic everywhere 
except perhaps at t = and t = oo. Because h{t) is single- valued it can only have power 
singularities (with integer exponent). Now 

1 /" „ , ^ , 1 /" / N dz , ^ 1 dz dz , , , 
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Figure 9: The eomplex t-plane. The contours corresponding to S and H divide the plane into 
sectors that correspond to different branches of the function g{z). The shaded regions correspond 
to gh{z) and (7r(5). (The interest of this picture hes in its quahtative features, but it was actuaUy 
computed from a numerical solution of the BAEs. The parameters are 0l — —0.46, (/ir = —0.653, 
riL = 152, nR = 186, and L = 200. These values correspond to 6 = 2ie^'^ °^'.) 



where the last integral is over some contour running from oo to 0, is finite. Since ^ remains finite 
and non-zero for t near or cxd, and 

dz _6{b~b*)t^ jt^ ift->0, 
dt~ {t^~ 1)2 ^ ift^oo, 

it follows that h{t) has at worst singularities at i = and at t = oo. Hence, the 1-form 

g{z) dz^h{t)%dt 
dt 

is nonsingular at i = and t — oo. 



5.4 Calculation of g{z) 



In the previous subsection it was shown that the 1-form (5.3) is nonsingular at t = and t = oo. 
The only singularities it can have are second order poles at the zeros ii, t2, . . . , of — 1. (These 
are the points in the i-plane corresponding to z = oo.) Therefore it can be written as 

The coefficients Vk and Sk are given by 

Tk = Rest^tj^ ^^^^~dl '^^ ^ R-<3S2=oo g{z)dz. 



and 



dz 

Sk = Rest=t^(i - tk)h{t)— dt = [[t - tk)z] Res^^oo z~^g{z)dz, 
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where the appropriate branch of g(z) is to be taken. 

The residues Res^^oo g{z)dz and Res^^oo z~^g{z)dz stiU have to be computed. From ( |60| ) and 
( |6l| ) one has 

Res£=oo5L(^)di = --^ / gni'n)dfi =: i?L, 
ReSi=oo 5R(^)d^ ^^TT^ I gUi) di =: i?R 

and 

Resf^oo z~'^g-L{z)dz = -1, 
ReSz=oo z^^9R{z)dz = -1. 

The residues for the other branches of g{z) follow by application of the monodromy operators. 
They are listed in Table ^. It follows from ( |50| ) that the integrals _Rl and i?R are real. 



Table 2: The poles and residues of g{z)dz. 



k 


tk 


5 


Resj^oo 9{z)dz 


ReSz=oo z ^g{z)dz 


1 


g7ri/3 


9L 




-1 


2 


_g-7ri/3 


-5R 


-Rr 


1 


3 


-1 


-5L - 5R 


— i?L — 


2 


4 




-5L 


— i?L 


1 


5 


g-iri/3 


9R 




-1 


6 


1 


5L + ffR 


-Rl + -Rr, 


-2 



Combining these results gives after some algebra that 



9i^) - E 



Sk 



t-tk [t- tkf 



= (1 - 2C)t + (1 - 2C*)t-i + C(< + + C*(<-i + t^) 



with 



6 2%/3Im6 



(64) 



We shall now argue in the generic case, b ^ 2i, that C = 0. From (|54| ) and (p5| ) the curves 2 
and H are described by Re[/L(2:) dz] = and Re[/R(z) dz] = 0, respectively, so the corresponding 
curves in the i-plane are both solutions of 



Re 



g{z) dz dz 
z dz dt 



dt 



0. 



(65) 



Note that z and g| are not single- valued functions of t, but the two branches of 



Idz 

z dz 



1 



differ only by a sign, which does not influence (^) . The two different solutions of (^5|) correspond- 
ing to S and H, respectively, meet at i = (and bX t = oo), so at these points the differential 
equation admits multiple solutions. When t — > 

9{z) dz dz ^ ^b-y_ ^ 



z dz dt 



b + b- 
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so this implies that C = 0. 

We shall now argue in the special case b — 2i that C — 0. When t 

/(^) dz = dt = 6 [cr' + (1 - cyt + o(<3)] dt 

(and similarly when t —> oo). The finiteness of the integral ( |63| ) (or its analogue for p^i) implies 
that C = 0. 

Now (|4|) becomes 

g(z)=t + i-i. (66) 

The functions 5l(-2) and 5r(2) are obtained by taking the appropriate branches ti^{z) and t-ii(z) 
of t{z). The branch iL(^) is determined by ti^{oo) = e'^'/^ and the fact that it has H as its only 
branch cut. Similarly tR(z) is determined by iR(oo) = c^'^^f-^ and the fact that it has S as its only 
branch cut. 



6 Calculation of physical quantities 

In Section ^ the relation was established between the canonical ensemble we are interested in 
and the semi-grand canonical ensemble that arises in the BA from Section |^. In Section |5| BA 
information was encoded in two complex functions satisfying a set of integral equations. These 
functions were then solved from those equations. In the present section the physical quantities 
occurring in Section ^ are extracted from the complex functions determined in Section ^j. 

6.1 Calculation of pL, Pr? 4>l, <Pr, and $ 



From ( |59| ) and (66) fhiz) and /r(z) are both given by 



= i±^^ (67) 



with different branches of t. It was claimed in Subsection 5.1 that the BA parameters pL, PR, '/'l 
and 0R and the semi-grand canonical free energy $ can be computed from the fmictions fh{z) 
and /r(z). They depend on the point b. The particle density pL was already computed in (p3): 



PL 



^Jj^{z)dz. (68) 



Because fi^{z) is analytic this integral does not depend on the precise shape of S, but on its 
homology only. 

Next (ji'L is calculated. Since fhiz) is known the function Fi^{z) is determined up to an inte- 
gration constant. The real part of this integration constant is fixed by ReFL(6L) = 0, see (|4^). 
From ( ^6| ) one has 

Re[FL(z) + FL(-z-i)] =2^L. 

It is now easy to compute that 

2 ^ ^ 



\e / Mz)dz. (69) 



From (P3| ) the free energy <I> equals — (Sl -I- Sr) with 

1 1 

El = lim y log |?7il, 

1 1 

Er= lim T^yXl^^SlC^I- 

i=l 
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Using ( [46D one calculates 

riR 

L 



Re [Fl (z) - log z] " = 2 - ^ log | , 



so 

EL^^Re^" (^/L(z)-i) dz. (70) 

In (^) and (^0|) the integral again only depends on the homology of the integration path, not on 
its precise shape. The real part of the integral even does not depend at all on the path chosen 
between the integration end points, but the imaginary part does. This is because the indefinite 



integral (46) is a sum of logarithms with real prefactors, and distinct branches of a logarithm differ 
by a multiple of 27ri, which is purely imaginary. 

Replacing all subscripts L in (|68|), ( |69| ) and ( [70| ) with R yields expressions of pR, and 
as integrals of functions involving /r(z). These integrals for pL, PR, 4'h, 0r, El and Sr are 
of the form Jydz, where the points {y,z) lie on an algebraic curve of genus 5. Therefore the 
indefinite integrals cannot be expressed in terms of "standard" functions. This does not prove 
that the definite integrals we need cannot be expressed in terms of standard functions, but it 
seems unlikely that it can be done. Of course they can be evaluated numerically. 

6.2 Calculation of |f and |f 



The Legendre transformation in Subsection 4.1 involves the derivatives and 7^-- Unfortu- 



nately we have not been able to compute $ as a function of pL, PR, and 0r for all values of these 



arguments. Instead we have in Subsection 6.1 computed these parameters and the free energy in 
the case that the curves S and H close, as a function of their common end point 5 = 6l = f'R- In 
order to still obtain the derivatives and we resort to a perturbation analysis. The details 
can be found in Appendix ^ here we only give some results. An infinitesimally small complex 
parameter C describes how far the curves open up. The thermodynamic parameters pL, Pr, 
(pn and the free energy $ then are functions of b and C. If all their first-order partial deriva- 
tives are known, and can be found by applying the standard coordinate transformation 

formula to the transformation between coordinates Re b, Im &, Re C and Im C on the one hand 
and PL, Pr, (/'l and (jj-n on the other hand. The derivatives with respect to Re 6 and Im6 can be 



obtained immediately from the integral expressions in Subsection 6.1. For the derivatives with 
respect to ReC and ImC the perturbation analysis is needed. It tells that to leading order in 
C the parameters pL, Pr, (f>L, <^r and the free energy $ = —(El -I- Er) are again given by the 



integrals (pq), (|69|), ( |70[ ) and their analogues involving f-R,{z), where f{z) is now given by 

z 

This yields integral expressions for their derivatives with respect to Re C and Im C. 

The expressions thus obtained for the partial derivatives of pL, Pr, 4>h, (^R and $ with respect 
to Re b, Im 6, Re C and Im C were evaluated numerically for some chosen value of 6, and from this 
and were calculated. These derivatives were also computed from numerical solutions of 
the BAEs for large system size L by numerical differentiation. The results from the two methods 
agree, which supports the perturbation analysis of Appendix ^ 

6.3 Configuration of S and H 

In the previous two sections several physical quantities have been expressed as integrals of functions 
involving /l(z) and /r(2:). These integrals depend on the parameter b and on the topology of the 
curves E and H, but not on their precise shape. If 6 ^ 2i there are two distinct points in the z-plane 
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corresponding to b, say bi and b2. The end points of S could be bi and 6J or 62 and 62, and the 
same holds for H. Therefore one can expect at least four different configurations for one and the 
same value of b. In order to determine what these four configurations actually are, we first guessed 
what they might look like. Then we chose some particular value of b (close to 2i) and for each of 
the four expected cases computed the value of the particle densities pL and and the phases 0l 
and 0R,. The BAEs were solved numerically for these parameter values, for large system size L. 
The resulting curves followed by the ^ and the rj display indeed the presupposed configurations. 
These curves are shown in Figure Note that without first guessing the configurations we would 
have had no way to find the values of the parameters pL, PR, 'P'l and so there would have been 
no BAEs to solve numerically. 





Figure 10: Four possible configurations of the curves S and H. The dashed curved are — 
and — H~^. In cases I and IV S and H have the same end points. In cases II and III S and — H^^ 
share end points, as do H and — S^^. 

The numerical results show that these four cases are related by the symmetries of the parameter 
space discussed in Subsection |4.2| . They are in a single orbit of the sub-group of order six generated 
by the horizontal translation (which is of order three) and the product of the reflection in a 
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horizontal line and the reflection in a vertical line (which is of order two) . For the remaining two 
members of this orbit we have not been able to numerically solve the BAEs. In these cases the 
particle densities are fairly high; we suppose that the curves S and H would cross or otherwise 
violate the condition that S and H only share their end points. 

Once this symmetry is known a (numerical) calculation of the physical quantities needs to be 
carried out only for one of the four cases I-IV. The values for the other three cases are then 
obtained at once by application of the symmetry transformations. 



6.4 Calculation of the sub-lattice densities and the entropy 



In Subsection 3.1 the physical quantities pL, PR, 0l, (/>r, and $ pertaining to the se mi-gr and 
canonical ensemble were calculated from the functions fhiz) and fK{z). In Subsection |6.2| the 
derivatives 



and 



were computed. Substitution of these results into formulas (36)-(|4l|) and 
( ^ ) from Subsection 4.1 gives the sub-lattice densities and the entropy, physical quantities for the 
canonical ensemble. This was performed numerically for a particular value of b. The results reveal 
that Po — P2 — Pa in cases I and IV and Pi = p^ = P5 in cases II and III. From the expressions 
(^)-(^l]) for the sub-lattice densities this is equivalent to 



5$ 



(2 + PL - 2pr) and 



9$ 



1 



d4>R 6 



(2 - 2pL + Pr) 



in cases I and IV and 



9$ 



2pr) 



and 



5* 

^(2pL - PR) 



dipK 6 



(72) 



(73) 



in cases II and III. 



One might hope to derive these expressions analytically from the results of 
Subsection |6.lH6.2|. We have not tried this because it would involve rather cumbersome relations 

Once the expressions (|7|) and (|7^) have been accepted 
becomes superfluous. Substituting them 



among integrals like (pSj), ( |69| ) and ([? 
the perturbation analysis approach from Subsection 

into (|3^)-(pl|) and (^) yields new expressions for the sub-lattice densities and the entropy. The 
expressions for the sub-lattice densities are polynomials in the particle densities pL and pR, the 
expression for the entropy also contains the phases 4>i, and (/)r and of course the free energy 

The cases I-IV correspond to different regions in the parameter space of sub-lattice densities, 
as given in Table ^. These four cases are defined for Re & > by Figure The mirror images 
(with respect to the imaginary axis) of the configurations in Figure 10 define cases I'-IV for 
Re 6 < 0. For example, the locus of 2 (H) for case I' is the mirror image of the locus of H (S) for 
case I. Table ^ also lists the regions in the parameter space of sub-lattice densities corresponding 
to the cases I'-lV. 



Table 3: The regions in parameter space of the sub-lattice densities for the cases I-IV and I'-IV'. 



I^lI 


I^rI 




ReS > 


Re6 < 


> 1 


> 1 


Po 


= P2 


= P4 


I 


pi > P5 > P3 


I' 


P5 > pi > P3 


> 1 


< 1 


pi 


= P3 


= P5 


II 


Po > P2 > Pi 


III' 


P4 > Po > P2 


< 1 


> 1 


Pi 


= P3 


= P5 


III 


P2 > Po > Pi 


II' 


Po > P4 > P2 


< 1 


< 1 


Po 


= P2 


= Pi 


IV 


P5> P3> Pi 


IV' 


Pi > P3 > P5 



6.5 Summary 



In the foregoing an exact solution of the trimer model was derived. Because the results are 
obtained in the course of a long derivation, we here provide a guide through them. The final result 
is the entropy as a function of six sub-lattice densities pi defined in Subsection 2.1. Com plete 
coverage of the lattice (|l|) and a further geometric constraint (||) (derived in Subsection 2^) leave 
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four independent parameters. A full analytic solution in the thermodynamic limit is obtained in 
a two-parameter subspace. 

The four-dimensional parameter space is described by new variables pL, PR, 0l and The 
free energy function of the ensemble with these parameters is denoted by $. The sub-lattice 



densities are expressed in pL, pr, ■§§- and in the entropy is given in terms of <i>, 



</>L, 0R, ^ and ^ in (^). 

The free energy is written as a sum, $ = — (Sl -f- Sr). In the solvable subspace the quantities 
PL, (/>L and El are expressed as contour integrals of a function /l(z) in (68)-(|7C|). Analogously 
the quantities pR, (/)r and Sr are integrals of a function fn{z)- The integration paths in the 
integral (|8|) for pL and its analogue for pR are contours S and H, respectively. These contours 
are symmetric under complex conjugation. Their end points in the upper half plane are denoted 
&L and 6r, respectively. These satisfy the equality 6l — = &r — &r^ = b. For each value of b 
this equation has two distinct solutions for 6l and for 6r , resulting in four configurations I-IV of 
S and H shown in Figure |lO|. The derivatives and are given by (^) in the cases I and 
IV and by (73) in the cases II and III. 

The functions f-L{z) and /r(2;) are different branches of a function J{z). The branch cuts of 
f-L{z) are H and — H^^, and /r(z) has branch cuts S and — S^^. In terms of a new variable t, 
defined by (^ ) and (62), the function /(z) is single- valued. It is given by (|67|), while the functions 
/l(2:) and /r(z) are reco vere d by selecting the appropriate branch ti^(z) and t^{z) of t, specified 
at the end of Subsection 5.4. 



7 Phase diagram 



In Subsections ^ and |2^ a linear and a quadratic constraint on the six sub-lattice densities were 
derived. In this section we first show that these constraints imply a breaking of the symmetry 
between certain sub-lattices. This symmetry breaking suggests that a phase transition takes place 
when the total density Pv = Pi + P3 + Ps of down trimers is increased from to 1. Next we 
compute the entropy as a function of pv from the exact solution of this model. From this entropy 
the phase diagram of the model in the parameter pv is obtained. It is also formulated in terms of 
the chemical potential of the down trimers instead of their density. 



7.1 Symmetry breaking 

The linear constraint (0) on the sub-lattice densities can be rewritten as 

po + P2 + P4 = 1 - Pv , (74) 
and from the quadratic constraint (||) one has 

P0P2 + P2P4 + P4P0 < ^Pv (75) 

with equality if and only if pi = ps = ps — ^Pv- If Pv is small (to be precise: smaller than 
2\/3 - 3), it follows from ^ and (|7|) that one of Po, p2 and p4 is larger than the other two, 
say Po > P2,P4. Thus the symmetry between the sub-lattices 0, 2 and 4 is broken. If there is no 
further symmetry breaking then p2 = p4 and pi ~ ps — ps, so 

Po > Pi = P3 = P5 > P2 = P4- (76) 

By he same token the symmetry between the sub-lattices 1, 3 and 5 is broken when pv is close 
to 1. When pv is increased from to 1 the following seems to be the simplest possible scenario. 
At Pv = sub-lattice is fully occupied and the other sub-lattices are empty. The six sub-lattice 
densities change continuously with pv, and ( |76| ) holds up to pv = ^. There the six sub-lattice 
densities are all equal to g. Then one of the odd sub- lattices, say 3, takes over, and 

P3 > Po = P2 = P4 > Pi = P5 
all the way to pv = 1 where all trimers sit on sub-lattice 3. 
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7.2 Entropy for pv 

In the previous subsection the occurrence was suggested of a phase transition when pv is increased 
from to 1. For the study of such a phase transition it would be helpful to know the entropy of the 
model as a function of pv = Pi +P3 + Ps- However, what we have computed thus far is the entropy 
as a function of all sub-lattice densities, but only for a two-dimensional subspace. Therefore for 
given pv the sub-lattice densities have to be determined for which the entropy is maximal. If we 
are fortunate these sub-lattice densities happen to lie in the two-dimensional solved subspace. 
For given pv < ^ the most symmetric possibility for the six sub-lattice densities is described 



by ([76|). Another possibility, 



P2 = Pi> Pi= pz= Pb> Po, 



(77) 



exists when 2\/3 — 3 < pv < 5- Because of symmetry ( |7^ ) and ( |77| ) are stationary points of the 
entropy. It is tempting to believe that ([76|), being the more general of the two most symmetric 
cases, corresponds to the maximum of the entropy. By numerically solving the BAEs the entropy 
of the model can be computed to high precision. Such calculations confirm that for pv < 5 the 
entropy takes its maximum at the symmetric case 
the solvable subspace. 



6[) of the sub-lattice densities, hence within 



As was seen in Subsection 3.4, for the solvable subspace one has po = P2 = Pi in cases I and IV 
and pi = P3 = P5 in cases II and III. Consider case II and take h on the imaginary axis between 
and 2i. The contours S and H then lie symmetric with respect to the imaginary axis, so pL = Pr, 
and hence p2 = Pi- Thus this is precisely the symmetric case (^). Therefore we have obtained 
the entropy as a function of pv for pv < 5. The entropy for pv > 5 follows immediately by the 
symmetry between up and down trimers. This entropy can also be obtained by considering case I 
and taking b above 2i on the imaginary axis. The resulting entropy is shown in Figure |ll|. When 
in case II h is not taken on the imaginary axis between and 2i, p2 7^ p4. Figure |l^ shows the 
entropy as a function of the asymmetry p2 — p4 at fixed pi = pa = ps along the line determined 
by the constraints and (|^). 




Figure 11: The entropy per trimer 5 as a function of the total density of down trimers pv = 
Pi + P3 + P5- It is obtained from the exact solution in the special case ( [tg] ) for pv < ^7 and 
similarly for pv > ^• 



For 6 = 2i all four cases I, II, III and IV coincide. The integrals in Subsection 6.1 then simplify. 
The sub-lattice densities are all equal to i, and the entropy per trimer is S'sym — log |\/3- 



7.3 Phase transition 

Consider a system with pv between and i. The energy is a convex function of pv, so the 
system is thermodynamically unstable. It would separate into a frozen phase with pv = and the 
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Figure 12: The entropy per trimer as a function of p2 — p4 at fixed pi = ps = P5- The end 
points of the lower curve (pv = 0.45) are determined by p2 = and by p4 — 0. The upper curve 
(pv — 0.48) terminates when = po and when p2 = Po- At the left end point of the middle curve 
(pv — 2\/3 — 3 w 0.4641) p2 = and p4 = po, at the right end point p4 = and p2 — po- 

symmetric phase with pv = 5- However an interface between these two phases is not possible in 
the model. Similarly a system with pv between ^ and 1 would demix into phases with pv = 
and Pv = 1, if coexistence between these phases were possible. 

Now give a chemical potential to the down trimers instead of imposing their density pv- 
The free energy 

F = -^vPv - 'S'(pv) 

takes its global minimum at 

{0 for Pv < -25sym 

i for -2S'syin < Pv < 25syin 
1 for 2S'sym < Pv 

Therefore the model is in a frozen phase for pv < — 2S'sym and for pv > 2S'sym and in the symmetric 
phase for — 2S'sym < Pv < 2S'sym- At pv — — 2S'sym and at pv = 2S'syni there is coexistence between 
a frozen and a symmetric phase. 

8 Conclusion 

We have introduced a new simple lattice model. It is a fluid of particles each occupying three 
sites of the triangular lattice. We distinguish six sub-lattices of adsorption sites for the trimers. 
Full occupancy and a resulting geometric constraint leave of the six sub-lattice densities four 
independent parameters. 

In the full four-dimensional parameter space the model is solvable by the Bethe Ansatz. In the 
thermodynamic limit the Bethe Ansatz equations can be reduced to two integral equations. In 
a two-dimensional subspace of the sub-lattice densities these integral equations can be solved by 
means of monodromy and analyticity properties of the functions involved. Within this subspace 
the entropy and the sub-lattice densities are given as integral expressions. 

The solution is very similar to that of the square-triangle random tiling model [^|l^. In 
both cases the solution is closely connected to the hexagonal domain wall structure of the model. 
Another solvable model with such a domain wall structure is the three-colouring model on the 
honeycomb lattice ||2^. In a configuration of this model the edges of the honeycomb lattice are 
coloured with three colours in such a way that the three edges meeting in each vertex have different 
colours. Alternatively this model can be formulated as the zero-temperature antiferromagnetic 
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three-state Potts model on the Kagome lattice [^l],^. We shall now briefly discuss the relation 
between these three models. 



The domain wall structure of the trimer model is depicted schematically in Figure 13, It 



contains two types of Y-joints but only one type of upside-down-Y-joints. In the square-triangle 
tiling there is only one type of Y-joints and one type of upside-down-Y-joints. In the honeycomb 
lattice three-colouring model on the other hand both the Y-joints and the upside-down-Y-joints 
come in two types. Hence these three models appear to be different. 



Figure 13: Schematic representation of the domain wall structure of the trimer model. The Y- 
joints of the domain walls come in two types. These are mirror images, featuring either a or 
a In contrast there is only one type of upside-down-Y-joints. 

The ^2^-' model is a vertex model on the square lattice, derived from an affine Lie alge- 
bra [^,p^. It satisfies the Yang-Baxter equation |^,|2^, so it can be solved by algebraic Bethe 
Ansatz [p7|. At a special value of the spectral parameter it is the three-colouring model on the 
honeycomb-lattice . For a suitable choice of the remaining parameters one of the two types of 
Y-joints and upside-down-Y-joints in the domain wall network is excluded. In this limit the model 
is just the square-triangle tiling. This mapping "explains" the solvability of the square-triangle 
tiling in terms of that of the A^'^ model [ p9| . 

In a similar fashion the square-triangle tiling can also be obtained from the trimer model. 
When the trimers on sub-lattice 4 (or 2) are excluded, one (or the other) type of Y-joint no longer 
occurs in the domain wall network. Again the square-triangle tiling results. The Bethe Ansatz 
for the trimer model remains valid in this special case. However, the substitutions (|2^) no longer 
makes sense when = (or W2 — 0), so the same is true of the analysis in the subsequent 
sections. 

Therefore the three models are connected in sense that both the trimer model and the A'^'' 
model contain the square-triangle random tiling as a singular limit. It would be interesting to 
know if the trimer model, like the square-triangle tiling, is a special case of some model satisfying 
the Yang-Baxter equation. 
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A Perturbation analysis 

In Subsection the quantities pi^, pn , (pi^, ipjx. and $ = — (Sl + Sr) were obtained as functions 
of 6 = 6l = For the computation of the sub-lattice densities and the entropy the derivatives 

and 



PlPr<Pl 



are also needed as functions of b. These cannot be calculated by differentiation of the $ already 
obtained, because variation of cjji^ (0r) at constant pL, PR and {(t>L) breaks the condition 

= ^R- Therefore in this appendix we infinitcsimally relax that condition, and compute pli PRi 

(pR and $ — + Sr) to leading order in the infinitesimal relaxation parameter. 

When the curves S and H do not close, the monodromy group is the full group SL(2,Z). 
Therefore gh{z) and 5R(i) are no longer single- valued functions of the variable t. Kalugin 
has provided a perturbation analysis for the analogous situation in the square-triangle random 
tiling model. It leans heavily on the understanding of the structure of the Riemann surface of the 
functions. Our approach does not require such knowledge and is more systematic. 

Although gh{z) and 5r(z) are no longer single- valued functions of the variable t, one can still 
perform the variable transformation (|6^). The end points 6l and S£ of S (6r and b'^ of H) then 
correspond to points cJl and rf^"^ (c^r and in the t-plane. The point b in (^) can be chosen 

such that |dL| = |c?f{ |; write 

dh = PlS and (Ir = PrS, 

with S real and positive and = |/3r,| = 1- The curves corresponding to S and H divide 
the annulus 6 < \t\ < into sectors, much as in Figure ^. We get a single-valued function 
h(t) = g{z) in this annulus instead of in the whole t-plane. Since it is analytic in the annulus it 
can be expanded as a Laurent series in t: 

oo 

g{z)^h{t)^ J2 AptP. (78) 

p— — C30 

Only powers with p = ±1 (mod 6) have the correct monodromy properties, so other powers 



cannot occur. From (50) one has 

hit*-')^g{z*)=9i^r = Ht)\ 

so the coefficients Ap satisfy 

A-p = a; . (79) 

We want to view the function g{z) given by ( fz^ ) as a perturbation of the function g{z) given 
by (|66|), where S is the small parameter. In our notation we have suppressed the dependence of 
the coefficients Ap on S, /3l, and /3r. 

The function g{z) satisfies the integral equation (|60|). We investigate how each of the terms tP 
from the Laurent series ( [ts] ) of g{z) behaves in this equation. In order to compute the integral we 
change from rj to t = t{fj) as integration variable. The resulting integrand is a rational function in 
T, which we decompose into partial fractions. Integration yields polynomial as well as logarithmic 
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terms; some care is required in choosing the branch of the logarithms. Finally we expand in powers 
of (5, obtaining 

2t:i J^rj- z ^ 27ri p - 6g 



g— — oo 



Y.^qit''-^)^''-'j (80) 

for each term in the Laurent series (|7^). Here t in the RHS corresponding to z in the LHS 
is in the sector containing ii, that is the sector where g{z) equals ghiz). Comparison with the 



integral equation (60) shows the following. The term in the RHS of (|80| ) exactly matches the 
term g-L{z) in the LHS of (|60|). The inhomogeneous term —t^ in the RHS of (|o|) corresponds 
to the inhomogeneous term 1 in the integral equation. The other terms in the RHS of ( |80| ) are 
"unwanted"; the powers of t they involve are multiples of 6. Because the Laurent series ( |78| ) 
satisfies the integral equation (|6^) the inhomogeneous terms —t^ from (|^ counterbalance the 
inhomogeneous term 1 of the integral equation, 

oc 

P— — 00 

(which means that (/l(oo) = 1), and the unwanted terms cancel, 

V -^SPA„ = forallg<0, (82) 
■^-^ p — bq 

P— — 00 

°° gP 

V -t^S-PA„ = forallq>0. (83) 
'^-^ p — oq 

P— — 00 

(Due to (^ the equations for q and for —q are equivalent.) The function g{z) also satisfies the 
integral equation (|6l|); this leads to another similar set of conditions on the coefficients Ap. 



The form of the equations ( |82| ) and ( p3D and their analogue from (61) suggests that for /3l and 
/?R fixed the coefficients Ap should be power series in S, 

Ap^A(^^+A^^^S + A(^^S' + ... (84) 

We would like to determine the coefficients A'^^ 

When t approaches the boundary of the annulus, \t\ —> S oi \t\ —>■ S^^ , the unperturbed function 
g{z) given by (|6^) becomes of the order S. It seems reasonable to assume that the terms Apt^ of 
the perturbed function g{z) do not grow faster than this, so the coefficients A^p''' with h < \p\ — 1 
must be zero. 

Consider (|8l| ) and its analogue from (|6l|). Substitution of the power series (|4|) yields, after 
rearrangement of the terms, 

oo / h+l \ 

E^M E for k = i,5. (85) 



The 5° part gives 



h=Q \p=-{h+l) 



tkA^"'' + t-^A^"l - 1 for fc = 1, 5. 



The unique solution of these equations reproduces the unperturbed function g{z) given by (|6q). 
For 1 < < 3 the S'' part of (p|) gives 



=0 for fc = l,5. 
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This implies that A[^^ and A''_'^1 are zero. The (5^ part of ( |85| ) gives 

tkA[^^ + t-'A^^l + 44'^ + t-'A^^l =0 for A; = 1, 5. 



These equations have two hnearly independent solutions, one of which satisfies (|79|). Substituting 
these results into (p3) and (fTSl) yields 



where we have written 



t + t-'^ + c (r^ ~t)+C* {t^ - t-^) + 0{S^) , 



(86) 



C 



and 



a'^^s' 



c*. 



Note that ( [Sq ) can be written in the form (64). We have used the equations (^ and (^ only to 
come up with the series expansion (|^. These equations would be needed if the coefficients A^'^ 
with h > A were to be determined. Knowledge of these coefficients would yield a solution to the 
integral equations (pO) and (61) also for a finite opening between &l and &R. Unfortunately we 



have not been able to calculate these coefficients, but fortunately they are not needed, because 
the present purpose is only to compute /Ol, /cr, 0l, '/'r and $ = — (Sl + Sr) only to leading order 
in 5. 

As our aim is to calculate the quantities pL, P R, (j^h^K and <&, we substitute ( [zsl) and (|8^ ) 



into the integral expressions from Subsection |6.1| . For (|6^) this gives, after transforming to t as 
integration variable: 



PL 



1 
27ri 



E 



E 



/3l<5 



1 



di 



^/|2T4d^ 



■At. 



(87) 



For each p and h we determine the order in 5 of the contribution. When t — )■ or t — > oo the 
integrand is proportional to t^'^^ and t^"^, respectively. Hence the integral is bounded, of order 
in 5 that is, for \p\ < 5; logarithmic in S for \p\ = 6; of order 6 — \p\ in S for |p| > 7. Note that the 

coefficients Ap''-' with |p| = 6 are zero. Let m denote the order in 6 of the integral. The order in S 
of the whole contribution is h + m: 



\p\ 


m 


h 


h + m 


1 











1 





4 


4 


1 





> 5 


> 5 


5 





4 


4 


5 





> 5 


> 5 


> 7 


6-bl 


> bl - 1 


> 5 



Therefore pL in (|83) has a 5" contribution from the unperturbed part in the RHS of (|8q), a 5"' 
contribution from the part involving C and C*, and contributions of higher order in 5 from the 
parts collected in the O(J^) term, so 



PL 



— / [t + t-^ + c - t) + C* (i^ - 1-^)] 



1 



dz 



/3L.5-1 



In the RHS the integration limits may be changed from f3i,S ^ and to oo and 0, as this makes 
a difference 0{S^). Transforming back to z as integration variable then gives 



PL 



/ [t + t-^+C (t-^ -t)+C* (t^ - t-^)] dz + 0(5'') , 



where ^^^^ denotes the unperturbed {5 — 0) contour. Hence pL is given to leading order in C ~ i^'* 
by (|68|), where now f{z) is given by (|7ll ) instead of (|67|), and integration is over the unperturbed 
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contour. Note that /3l and /3r do not occur in this expression. Similar arguments show that fuUy 
analogous results hold for pR, 0l, 0r, and Sr: up to 0{S^) they are given by the integrals 
(|68|), ( |69| ) and (^0|) or their R-analogues, with f{z) given by (|7l|). Therefore we have now obtained 
these quantities to leading order, namely S^, in the parameter S that describes the infinitesimally 
small opening 6l — between the end points of S and H. 
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